#!/usr/bin/perl
# 2013/08/19 Convert association file to BED format to draw in IGB

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 8 or $ARGV[0] eq '--help'){
	die "Usage : $0 -m [map file] -o [bed format] -n [name] -a [association file]\n";
}

my %opt;
getopts("m:o:n:a:", \%opt);
my $FILE_map  = $opt{m};
my $FILE_aso = $opt{a};
my $FILE_bed = $opt{o};
my $BED_name = $opt{n};


#browser position chr1:1000000-1005000
#track name="sampleBed" description="sample of ChIA-PET drawoing" itemRgb="On"
#chr1 1000100 1000900 none 0 + 1000100 1000900 255,0,0 2 25,25 0,775
#chr1 1000100 1000910 none 0 + 1000100 1000910 255,0,0 2 25,25 0,785
#chr1 1000100 1000910 none 0 + 1000100 1000910 255,0,0 2 25,25 0,785
#chr1 1000100 1000910 none 0 + 1000100 1000910 255,0,0 2 25,25 0,785
#chr1 1000100 1000910 none 0 + 1000100 1000910 255,0,0 2 25,25 0,785


#---------------------------------------
# read association
#---------------------------------------
my %association;
my $fh_aso = IO::File->new($FILE_aso) or die "cannot open $FILE_aso: $!";
while($_ = $fh_aso->getline()){
	s/\r?\n//;
	my ($peak1, $loc1, $peak2, $loc2, $interaction, $pval, $qval) = split /\t/;
#	if($qval eq 'NA'){
#		next;
#	}
#	unless($qval < 0.05){
#		next;
#	}

	unless($pval < 0.05){
		next;
	}

	my ($chr1, $pos1) = split /:/, $loc1;
	my ($start1, $end1) = split /\-/, $pos1;
	my ($chr2, $pos2) = split /:/, $loc2;
	my ($start2, $end2) = split /\-/, $pos2;

	if($chr1 ne $chr2){
		next;
	}
	if(abs($start1 - $start2) > 100000){
		next;
	}
	$association{$peak1}{$peak2} = 1;
}
$fh_aso->close();



my %data;
# id chr1 loc1 direction1 quality1 hind1 hind_loc1 chr2 loc2 direction2 quality2 hind2 hind_loc2 type(unique or repeat) number(alignable locations) locations
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $peakName1, $peakLoc1, $chr2, $loc2, $direction2, $mapQ2, $peakName2, $peakLoc2, $type, $num, $locations) = split /\t/;

	if($chr1 ne $chr2){
		next;
	}

	if(abs($loc1 - $loc2) > 50000){
		next;
	}

	if(abs($loc1 - $loc2) < 1000){
		next;
	}

	unless((exists $association{$peakName1}{$peakName2}) or (exists $association{$peakName2}{$peakName1})){
		next;
	}

	if($loc2 < $loc1){
		($loc1, $loc2) = ($loc2, $loc1);
	}

	$data{$chr1}{$loc1}{$loc2} = sprintf("$chr1 $loc1 $loc2 none 0 + $loc1 $loc2 255,0,0 2 25,25 0,%d", $loc2-$loc1-25);
}
$fh_map->close();


my $fh_out = IO::File->new($FILE_bed, 'w') or die "cannot write $FILE_bed: $!";
my $headder =<<"EOL";
track name="$BED_name" description="association of $BED_name" itemRgb="On"
EOL
$fh_out->print("$headder");
foreach my $chr(qw(I II III)){
	foreach my $i(sort {$a <=> $b} keys %{$data{$chr}}){
		foreach my $j(sort {$a <=> $b} keys %{$data{$chr}{$i}}){
			$fh_out->print("$data{$chr}{$i}{$j}\n");
		}
	}
}

$fh_out->close();